Effect of External Normal and Parallel Electric Fields 
on 180° Ferroelectric Domain Walls in PbTiOs 
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We impose uniform electric fields both parallel and normal to 180° ferroelectric domain walls in 
PbTiOs and obtain the equilibrium structures using the method of anharmonic lattice statics. In 
addition to Ti-centered and Pb-centered perfect domain walls, we also consider Ti-centered domain 
walls with oxygen vacancies. We observe that electric field can increase the thickness of the domain 
wall considerably. We also observe that increasing the magnitude of electric field we reach a critical 
electric field E"^; for E > E'^ there is no local equilibrium configuration. Therefore, E'^ can be 
considered as an estimate of the threshold field Eh for domain wall motion. Our numerical results 
show that Oxygen vacancies decrease the value of E'^. As the defective domain walls are thicker than 
perfect walls, this result is in agreement with the recent experimental observations and continuum 
calculations that show thicker domain walls have lower threshold fields. 



I. INTRODUCTION 

Ferroelectric materials have been used in many impor- 
tant applications such as high strain actuators, electro- 
optical systems, non- volatile and high density memories, 
etcjii^. The properties of domain vifalls in ferroelectric 
materials including their structure, thickness, and mo- 
bility are important parameters as they determine the 
performance of devices that use these materials^. 

Theoretical calculations have predicted that ferroelec- 
tric domain walls are atomically sharp and their thick- 
ness is about a few angstroms^"—. However, experimental 
measurements show the existence of domain walls with 
thicknesses of a few micrometers^iS. It has been observed 
that such broadening of domain walls is due to the pres- 
ence of extrinsic defects, charged walls, and surfaces^*'. 
Shilo et al.^^ used atomic force microscopy to measure the 
surface profile close to emerging domain walls in PbTiOa 
and then fitted it to the soliton-type solution of GLD the- 
ory. They measured wall widths of 1.5nm and 4nm and 
observed a wide scatter in wall widths. They suggested 
that the presence of point defects is responsible for such 
wide variations. Lee et al^ proposed a continuum model 
to investigate this proposal and reproduced the experi- 
mentally observed range of wall widths with their model. 
They mentioned that the interaction between the order 
parameter and point defects and interaction of point de- 
fects with each other are two important interactions that 
should be considered properly in such modelings. Jia et 
al^ investigated the cation-oxygen dipoles near 180° do- 
main walls in PbZro.2Tio.8O3 thin films. They measured 
the width and dipole distortion across domain walls us- 
ing the negative spherical-aberration imaging technique 
in an aberration-corrected transmission electron micro- 
scope and observed a large difference in atomic details 
between charged and uncharged domain walls. 

External electric field can cause the motion of ferroelec- 
tric domain walls if the magnitude of the field reaches the 
threshold field Eh for wall motion, i.e., the field at which 
a domain wall begins to move after overcoming the in- 
trinsic Peierls friction of the ferroelectric lattice-". It was 



observed that threshold fields that are predicted via ther- 
modynamic calculations are usually much greater than 
the experimental values. For example, Bandyopadhyay 
and Ray^'^ predicted an upper limit for Eh of LiNbOa 
to be 30000fcy/cTO but experimental observations show 
that the threshold field for wall motion can be less than 
I'akV/cra. Choudhury et ai— suggested that the reason 
for such large differences between theoretical and exper- 
imental values of Eh is broadening of the domain walls. 
Using microscopic phase-field modeling, they show that 
the threshold field for moving an antiparallel ferroelectric 
domain wall dramatically drops by two or three orders of 
magnitude if the wall was diffused by only about 1 — 2nm. 
Su and Landis^^ developed a continuum thermodynamics 
framework to model the evolution of ferroelectric domain 
structures and investigated the fields near 90° and 180° 
domain walls and the electromechanical pining strength 
of an array of line charges on these domain walls. 



In this work, we investigate the effect of external elec- 
tric field (E) on the perfect and defective 180° domain 
walls in PbTiOa using the method of anharmonic lattice 
statics. We consider both Pb-centered and Ti-centered 
perfect domain walls and also defective domain walls with 
oxygen vacancies. In agreement with experimental re- 
sults, our calculations show that such defective domain 
walls are thicker than perfect wallai^. By increasing E 
we reach a critical value E'^ such that for E > E'^ the 
lattice statics iterations do not converge. Therefore, this 
critical value can be considered as a lower bound for the 
threshold field for wall motion. 



This paper is organized as follows. In ilTTl we explain 
the geometry of the perfect and defective domain walls 
that we use throughout this work. In mill we describe 
the method of analysis used in our calculations. Our 
numerical results are presented in ijlVI The paper ends 
with some concluding remarks in ^|Vl 
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FIG. 1: (a) The relaxed configuration of the unit cell of PbTiOs. a and c are the tetragonal lattice parameters. Note that 
01, 02, and 03 refer to oxygen atoms located on (001), (100), and (OlO)-planes, respectively. (5 denotes the y-displacements of 
the atoms from their centerosymmetric positions and arrows near each atom denote the direction of these displacements, (b) 
The geometry of a perfect Ti-centered 180° domain wall, (c) The geometry of an Ol-defective 180° domain wall. Note that 
Pb-centered domain walls with oxygen vacancies are not stable. 



II. FERROELECTRIC DOMAIN WALLS 

Due to the relative displacements between the center 
of the positive and negative charges, each unit cell of 
a ferroelectric crystal has a net polarization below its 
Curie temperature. FiglIJa) shows the relaxed unit cell 
of tetragonal PbTiOa. In this work, we consider 180° do- 
main walls in PbTiOa parallel to a (lOO)-plane. These do- 
main walls are two dimensional defects across which the 
direction of the polarization vector switches. There are 
two types of perfect 180° domain walls in PbTiOa: Pb- 
centered and Ti-centered domain walls. Figlljb) shows 
the geometry of a Ti-centered domain wall. 

In addition to perfect domain walls, we also consider 
180° domain walls with oxygen vacancies. It is known 
that oxygen vacancies tend to move toward domain walls 
and pin themi^"— . Therefore, we study domain walls 
with oxygen vacancies sitting on them. In order to be 
able to obtain a solution, we need to consider periodi- 
cally arranged vacancies on the domain walls. Although 
in reality oxygen vacancies have lower densities, our re- 
sults with the current assumption can still provide im- 
portant insights on the effect of oxygen vacancies on 
180° domain walls. Depending on which oxygen in the 
PbTiOa unit cell sits on the domain wall, there would be 
three types of defective domain walls: (i) 02-defective, 



(ii) Ol-defective, and (iii) 03-defective. FiglD^c) shows 
an Ol-defective domain wall. Note that 01- and 03- 
defective domain walls are Ti-centered while 02-defective 
domain wall is Pb-centered. It has been observed that 
02-defective domain walls are not stabl e^^i^^ , i.e., the 
lattice statics iterations do not converge. Thus, we con- 
sider 01- and 03-defective domain walls in the following. 

Let x, y, and z denote coordinates along the (100), 
(010), and (OOl)-directions, respectively. We assume 
a ID symmetry reduction, which means that all the 
atoms with the same x-coordinates have the same dis- 
placements. Therefore, we partition the 3D lattice L as 
^ — U/ Uqgz 'C/q,, where and TL are 2D equivalence 
classes parallel to the (100) plane and the set of integers, 
respectively, j = J/3 is the atom in the fith equivalence 
class of the Jth sublattice. See^^— for more details on 
the symmetry reduction. 



III. METHOD OF CALCULATION 

We apply a uniform electric field on 180° domain walls 
and obtain the equilibrium structure using the method 
of anharmonic lattice statics^^ . We use a shell potential 
for PbTiOa -1 for modeling the atomic interactions. Each 
ion is represented by a core and a massless shell in this 
potential. Let £ denote the collection of cores and shells. 
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I € £ denotes a core or a shell in £, and {x*} .^^ repre- where C = (H* 
sents the current position of cores and shells. In this shell 
potential, three different energies are assumed to exist 
due to the interactions of cores and shells: i^gjjort' ^long' 

and ^core-shell- ^short {{'^'}^Gc) denotes the energy of 
short range interactions, which are assumed to be only 
between Pb-0, Ti-0, and 0-0 shells. The short range in- 
teractions are described by the Rydberg potential of the 
form {A + Br) exp(— r/C), where A, B, and C are poten- 
tial parameters and r is the distance between interacting 

elements, ^long lier) denotes the Coulombic in- 
teractions between the core and shell of each ion with 
the cores and shells of all the other ions. For calculat- 
ing the classical Coulombic energy and force, we use the 

damped Wolf method^^. Finally, ^^core-shell 
represents the interaction of core and shell of an atom 
and is assumed to be an anharmonic spring of the form 
{l/2)k2r'^ + (l/24)fc4r^, where ^2 and fc4 are constants. 
The total static energy is written as 



i+l 



V£\ and 



C • A 



aT . C* • a' 



(5) 



^ ^ec) = Wt {{^'}rec) +%ng ({^'}.e£ 

+ '^core-sheU ({^'}»e£) ■ (1) 



Note that all the calculations are done for absolute zero 
temperature. At this temperature PbTiOs has a tetrag- 
onal unit cell with lattice parameters a — 3.843 A and 
c = 1.08a^. 

Assume that a uniform electric field E — {E^, Ey,Ez) 
is applied to a collection of atoms. Then for the relaxed 
configuration B = {x'}jg£ C R-^, we have 



Vie/;, 



(2) 



where qi denotes the charge of the ith particle (core or 
shell). To obtain the solution of the above problem, we 
utilize the Newton method. Having a configuration B'' 
the next configuration S'^+i is calculated from the cur- 
rent configuration B'' as: 6*^+1 = B'' + 5^ , where 



5'^ = -h-mbM -VS B" 



(3) 



with H denoting the Hessian matrix. The calcula- 
tion of the Hessian becomes inefficient as the size of 
the problem increases and hence we use the quasi- 
Newton method. This method uses the Broyden- 
Fletcher-Goldfarb-Shanno (BFGS) algorithm to approx- 
imate the inverse of the Hessian^'^ instead of the direct 
calculation of the Hessian at each iteration. We start 
from a positive-definite matrix and use the BFGS algo- 
rithm to update the Hessian at each iteration as follows: 
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(4) 



Calculating C*+^, one then should use C'+^ instead of 
to update the current configuration. If C'+^ is 
a poor approximation, then one may need to perform 
a linear search to refine S^-'+i before starting the next 
iteration^'^. 

In the presence of oxygen vacancies on the domain wall, 
one needs to consider charge redistribution between some 
ions. To model an oxygen vacancy using a shell poten- 
tial, we remove the core and shell of the oxygen atom 
and because we assume a charge neutral oxygen vacancy, 
there will be a charge redistribution in the neighboring 
shells^-. It is known that charge redistribution is highly 
localized and hence in our calculations we equally dis- 
tribute the charge AQ = Qs + Qc, where Qs and Qc are 
oxygen shell and core charges, between the (fourteen) 
first nearest neighbors of each oxygen vacancy. 

To obtain the equilibrium configuration under an ex- 
ternal electric field we need to start from an appropri- 
ate initial configuration. This initial configuration for 
perfect and defective domain walls is the equilibrium 
configuration of these domain walls under zero electric 
field (see^'^*' for discussions on how to calculate these 
configurations). As we mentioned earlier, we assume a 
ID symmetry reduction for the lattice and hence as is 
shown in Fig[2l our computational box (CB) consists of 
a row of unit cells perpendicular to the domain wall. In 
this figure, the shaded region is the computational box. 
Note that because in general there is no symmetry in 
the problem, we need to relax all the atoms inside the 
CB. For removing the rigid body translation freedom of 
the atoms, one should fix the core of an atom and relax 
the other atoms. We fix Pb-core (Ti-core) of an atom 
located on the domain wall in Pb-centered (Ti-centered) 
domain walls. Thus, if there are M unit cells in the CB, 
we would have 30M — 3 variables in our calculations. We 
should mention that to investigate the effect of the size 
of CB in the domain wall plane, we consider CBs with 
one, four, and sixteen unit cells in the domain wall plane 
and therefore the number of the unit cells in CB in each 
case is M, 4M, and 16M, respectively. We observe that 
the final relaxed structure does not depend on the size 
of CB in the domain wall plane. This suggests that the 
symmetry reduction that we use in our calculations is a 
reasonable assumption for this problem. 

Note that we consider a finite number of unit cells in 
the CB and do not assume any periodicity condition in 
our calculations. This means that we need to impose 
some proper boundary conditions to take into account 
the effect of the atoms located outside of CB. To this 
end, we rigidly move the unit cells outside of CB with 
displacements equal to those of the first or last unit cell 
of the CB (the unit cell on the boundary of the CB that is 
closer to the unit cell outside of the CB). This is a natural 
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FIG. 2: Ti-cores under external electric field in a Ti-centered 180° domain wall. Ex and Ey axe the normal and parallel 
electric fields, respectively. The shaded region denotes the region that is relaxed in each step. Note that M is the size of the 
computational box (CB) normal to the domain wall. We consider different CBs with one, four, and sixteen unit cells in the 
domain wall plane. 



boundary condition as we expect the bulk configuration 
far from the domain wall. 

In our calculations we use M = 20 as larger values 
for M do not affect the results. Imposing an external 
electric field should be done step by step, i.e., one first 
needs to obtain the configuration for E = AEi from the 
initial configuration and then use this configuration to 
obtain the equilibrium configuration for E = AEi + AE2 
and so on. We use the average step size of 20kV/cm for 
electric field. Using this step size and force tolerance of 
O.OOSeFA , our solutions converge after about 30 to 40 
iterations. 



IV. NUMERICAL RESULTS 

In this section we present our numerical results for per- 
fect and defective domain walls. Note that as the coordi- 
nates of cores and shells are close to each other, we only 
report the results for cores. Also as we mentioned ear- 
lier, X, y, and z are coordinates along the (100), (010), 
and (OOl)-directions, respectively. 

Perfect domain walls: We plot the y-coordinates of 
Ti-cores under external electric field normal to the Ti- 
centered domain wall. Ex, in Figl3l[a). As expected, we 
see that increasing the electric field, the atomic struc- 
ture loses its symmetry. We observe that there exists an 
upper bound for E^-, i.e., there exists a critical electric 
field E^ such that for E^ > E^ there is no local equilib- 
rium structure. The critical value of the normal electric 
field is about E!^ = UOOkV/cm. The thickness of the 
domain wall slightly increases as the normal electric field 
increases. Note that domain wall thickness cannot be de- 
fined uniquely very much like boundary layer thickness 



in fluid mechanics. Here, domain wall thickness is by 
definition the region that is affected by the domain wall, 
i.e. those layers of atoms that are distorted. One can 
use definitions like the 99%-thickness in fiuid mechanics 
and define the domain wall thickness as the length of 
the region that has 99% of the far field rigid translation 
displacement. What is important here is that no mat- 
ter what definition is chosen, domain wall "thickness" 
increases as the normal electric field increases. For a 
Ti-centered domain wall, the domain wall thickness in- 
creases from 3 atomic spacings (Inm) to about 5 atomic 
spacings (1.5nm) for E^ = E^. 

FiglSlJb) depicts the y-coordinates of Ti-cores under 
an external electric field Ey parallel to a Ti-centered do- 
main wall. It is observed that such electric fields do 
not alter the domain wall thickness. Note that simi- 
lar to the atomic structure for normal fields, the atomic 
structure under parallel fields loses its symmetry as well. 
The critical value of the parallel electric field is about 
Ey — 5900kV/ cm, which is 4 times larger than that of 
the normal electric field. 

FigEJc) shows the y-coordinates of Pb-cores of a Pb- 
centered domain wall under normal electric field E^ . We 
observe that the critical electric field is about E^ = 
6300kV/cm, which is about 4.5 times greater than the 
critical normal electric field of Ti-centered walls. Also it 
is observed that domain wall thickness increases to about 
11 atomic spacings (Anm) under critical normal electric 
field. The y-coordinates of Pb-cores of a Pb-centered 
domain wall under parallel electric field Ey are shown in 
FigOl^d). Similar to perfect Ti-centered domain walls, we 
observe that parallel electric fields do not affect the do- 
main wall thickness. The critical parallel electric field is 
about Ey = 6500kV/ cm. For Pb-centered domain walls 
we see that unlike Ti-centered domain walls, the critical 
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FIG. 3: The y-coordinates of cores under external electric field: Ti cores in a perfect Ti-centered domain wall under (a) Ex 
and (b) Ey\ Pb cores in a perfect Pb-centered domain wall under (c) Ex and (d) Ey. 



normal electric field is close to the critical parallel electric 
field. 

Fig|4] depicts the polarization profiles normal and par- 
allel to the domain walls. For calculation of the cell- 
by-cell polarization, we follow Meyer and Vanderbilf*. 
We plot P = {Px,Py) = P/IPfcl, where P is the po- 
larization and jPbl — is the norm of the 
bulk polarization^^. FigUJ^a) shows and Py for a Ti- 
centered domain walls under zero and critical electric 
fields. In agreement with Lee et al^ and Angoshtari 
and Yavarii, it is observed that (100) Ti-centered domain 
walls have a mixed Ising-Neel character, i.e., polarization 
rotates normal to the (lOO)-plane near the domain wall. 
For E — the maximum normal component of the polar- 
ization is about 2% of the bulk polarization. For E = E^, 
as can be expected, normal electric field causes the pos- 
itive and negative charges to have normal displacements 
that create a polarization in the x-direction. This nor- 
mal component of the polarization (Px) reaches to about 
13.5% of the bulk polarization at E^, but we observe that 
normal electric field E!: does not have a remarkable ef- 



fect on the parallel component of polarization, Py. On 
the other hand, we observe that under E = Ey, P^ does 
not change considerably but Py has an unsymmetric pro- 
file with the maximum value of about 105% of the bulk 
polarization. 

FigHJb) presents similar results for Pb-centered do- 
main walls. Similar to Ti-centered domain walls, we ob- 
serve that Pb-centered domain walls have a mixed Ising- 
Neel characteri2i^ with P^ about 2% of the bulk polar- 
ization for zero electric field. For E — E^, P^ reaches 
to about 38% of the bulk polarization. Also we observe 
that E^ has more impact on Py compared to Ti-centered 
walls. Finally, it is observed that similar to Ti-centered 
domain walls, Ey does not have a significant effect on 
Px but makes Py unsymmetric with maximum value of 
about 107% of the bulk polarization. 

Defective domain walls: In this part we report the 
structure of defective domain walls under normal and 
parallel external electric fields. Because the results for 
01- and 03-defective domain walls are similar, we only 
present the results for 01-defective walls, which are Ti- 
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FIG. 4: The polarization profiles P = {Px,Py) of domain walls under zero, normal critical field (E^), and parallel critical field 
(Ey) for (a) Ti-centered, and (b) Pb-centered domain walls. 
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FIG. 5: The y-coordinates of Ti cores in an 01-defective domain wall under (a) normal (E^), and (b) parallel (-Ey) external 
electric fields. 



centered. Note that as we mentioned earlier, 02-defective 
domain walls, which are Pb-centered, are not stable. Our 
calculations show that they are not stable even under 
external electric fields. We had earlier shown that they 
are not stable under strain as welU^. 

Figinija) depicts the y-coordinates of Ti-cores in an 01- 
defective domain wall under normal external electric field 
E^. The critical normal field is about = 380kV/cm. 
It is observed that domain wall thickness increases up to 
about 16 atomic spacings {6nm) under the critical normal 
electric field. Comparing 01-defective atomic structure 
with the structure of perfect Ti-centered domain wall un- 
der normal field (FiglSja)), we observe that oxygen va- 
cancies increase the thickness of the domain wall con- 
siderably. Also it is observed that critical normal electric 
field of defective domain walls is smaller than that of per- 



fect Ti-centered wall. FiglSjb) shows the y-coordinates of 
Ti-cores in an 01-defective domain walls under parallel 
electric field Ey. The value of the critical field is about 
Ey — 51Q0kV/cm. Here we observe a major difference 
between the atomic structures of perfect and defective 
domain walls; unlike perfect domain walls, parallel elec- 
tric fields increase the thickness of defective domain walls 
up to about 13 atomic spacings (5nm) under the criti- 
cal parallel electric field. Also similar to normal electric 
fields, we observe that critical electric field of defective 
domain walls is smaller than that of perfect domain walls. 

Defective domain walls are thicker than perfect domain 
walls. The observation that the defective domain walls 
have smaller critical electric fields is in agreement with 
the experimental observations of Choudhury et alr^. 
They observed that the threshold field for domain wall 
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FIG. 6: The polarization profiles of 01-defective domain walls 
under (a) zero, (b) normal critical field (E^), and (c) parallel 
critical field (Ey). (d) Components of the polarization vector 
P = {Px,Py) of 01-defective domain walls under zero, -BJ, 
and El. 



motion exponentially decreases as the domain wall width 
increases. 

Figl6] shows the polarization profiles for 01-defective 
domain walls. It is observed that similar to perfect do- 
main walls, defective domain walls have an Ising-Neel 
character with of about 2.5% of the bulk polarization 
for zero electrical field. For E — E^, P^ reaches to about 
55% of the bulk polarization, which is greater than the 
corresponding values for perfect domain walls, and Py 
shows more Ising-type character. As we mentioned ear- 



lier, for E = Ey we observe a difference between perfect 
and defective domain walls; unlike perfect walls, parallel 
electric fields have considerable effects on P^: it reaches 
to about 55% of the bulk polarization under Ey. Similar 
to perfect domain walls, Py has an unsymmetric distribu- 
tion and reaches to about 106% of the bulk polarization. 



V. CONCLUDING REMARKS 

In this work we obtained the atomic structure of per- 
fect and defective 180° domain walls in PbTiOs under 
both parallel and normal external electric fields using the 
method of anharmonic lattice statics. We observe that 
electric field can increase the thickness of a domain wall 
considerably (up to 5 times thicker than domain walls 
under no external electric field). This can be one rea- 
son for the wide scatter of the domain wall thicknesses 
observed in experimental measurements. In agreement 
with previous worksii^i^, we observe that oxygen vacan- 
cies can increase the thickness of the domain walls. We 
also observe that by increasing the external electric field 
we reach a critical electric field E'^. For E > E'^ there 
is no local equilibrium configuration and hence E'^ can 
be considered as an estimate of the threshold field for 
the domain wall motion. We observe that defective do- 
main walls, which are thicker than perfect domain walls, 
have smaller critical fields. This is in agreement with the 
experimental observations that show the threshold field 
decreases as the domain wall thickness increasesiSi. 

In practice, it has been observed that domain walls 
move or break down under electric fields in the order of 
a few kV/cmi^:^^, which are considerably smaller than 
the high-fields that we consider here. We do not consider 
break down of the domain walls in our model. Also as was 
mentioned earlier, the high density of oxygen vacancies 
that we assume is unrealistic. In practice, steps and other 
complex defects on domain walls can increase the thick- 
ness of the domain walls considerably^i2i2i. Therefore, 
as the threshold fields for domain walls decrease expo- 
nentially with the increase of the domain wall widtbi^, 
one can obtain better estimates for the critical electric 
fields with more realistic models for defects on domain 
walls. Also as suggested by Roy et alH^ electric fields 
change the potential parameters. In this paper our aim 
is to demonstrate that even with our simple model, one 
can show that the threshold field has an inverse relation 
with the domain wall thickness. 



Acknowledgments 

We benefited from a discussion with CM. Landis. We 
thank an anonymous reviewer whose comments helped 
us improve the paper. 



8 



* Electronic address: 'arash.yavari@ce. gatech.edu| 
^ J. F. Scott, Science, 315 , 954 (2007). 

^ S. V. Kalinin, A. N. Morozovska, L. Q. Chen and B. J. 

Rodriguez, Rep. Prog. Phys., 73 , 056502 (2010). 
^ C-L. Jia, S-B. Mi, K. Urban, I. Vrejoiu, M. Alexe, D. Hesse, 

Nature Mater., 7, 57 (2008). 

* B. Meyer and D. Vanderbilt Phys. Rev. B, 65, 104111 
(2002). 

^ J. Padilla, W. Zhong and D. Vanderbilt Phys. Rev. B, 53, 
R5969 (1996). 

® A. Yavari, M. Ortiz and K. Bhattacharya, Philos. Mag., 
87 , 3997 (2007). 

A. Angoshtari and A. Yavari, EPL, 90, 27007 (2010). 
® M. Iwata, K. Katsuraya, I. Suzuki, M. Maeda, N. Yasuda 

and Y. K. Ishibashi, Jpn. J. Appl. Phys., 42, 6201 (2003). 
® P. Lehnen, J. Dec and W. Kleemann, J. Phys. D: Appl. 

Phys., 33, 1932 (2000). 
^° S. Choudhury et al., J. Appl. Phys., 104, 084107 (2008). 

D. Shilo, G. Ravichandran and K. Bhattacharya, Nature 

Mater., 3, 453 (2004). 

W. T. Lee, E. K. H. Salje and U. Bismayer, Phys. Rev. B, 
72, 104116 (2005). 

A. K. Bandyopadhyay and P. C. Ray, J. Appl. Phys., 95, 
226 (2004). 

Y. Su and C. M. Landis, J. Mech. Phys. Solids, 55, 280 
(2007). 

L. X. He and D. Vanderbilt, Phys. Rev. B, 68, 134103 



(2001). 

M. Calleja, M. T. Dove and E. K. H. Salje, J. Phys.: Con- 
dens. Matter, 15, 2301 (2003). 

L. Goncalves-Ferreira, S. A. T. Redfern, E. Artacho, E. 
Salje and W. T. Lee, Phys. Rev. B, 81, 024109 (2010). 
A. Angoshtari and A. Yavari, Comput. Mater. Sci., 48, 
258 (2010). 

A. Yavari, M. Ortiz and K. Bhattacharya, J. Elasticity, 86 
, 41 (2007). 

A. Yavari and A. Angoshtari, Inter. J. Solids Struct., 47, 
1807 (2010). 

A. Asthagiri, Z. Wu, N. Choudhury and R. E. Cohen, Fer- 
roelec, 333 , 69 (2006). 
22 D. P. Wolf, P. Keblinski, S. R. Phillpot and J. Eggebrecht, 
J. Chem. Phys., 110, 8254 (1999). 

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. 
Flannery, Numerical recipes: the art of scientific comput- 
ing (Cambridge University Press, Cambridge, 1989). 
2* A. Angoshtari and A. Yavari, J. Appl. Phys., 108, 084112 
(2010). 

25 D. Lee, R. K. Behera, P. Wu, H. Xu, Y. L. Li, S. B. Sinnott, 
S. R. Phillpot, L. Q. Chen and V. Gopalan, Phys. Rev. B, 
80, 060102(R) (2009). 

A. Roy, M. Stengel, and D. Vanderbilt, Phys. Rev. B, 81, 
014102 (2010). 



